The main data used in this tutorial and in the lecture are about the geolocalisation of french restaurants in Paris and in a department called Haute-Garonne. We use two different sources:
SIRENE has the advantages of being rigorous and exhaustive on the French territory.
OSM has many benefits, ensuring transparent data provenance and ownership, enabling real-time evolution of the database and, by allowing anyone to contribute, encouraging democratic decision making and citizen science.
sf::st_read function also work? Why?Use the readRDS function.
It would not work because ‘iris_31’ is not a shapefile but a file already R formatted. Simply load it with the readRDS function.
plot(iris_31). What do you notice ?We notice that R performs 3 graphs: one graph per variable in the sf object.
sf::st_geometry function? What solution do you propose then?sf::st_geometry makes it possible to isolate the information contained in the ‘geometry’ column of the sf object. Using it, we put aside other variables (here CODE_IRIS, P14_POP and INSEE_COM).
Test the Azimuthal Equidistant projection with crs="+proj=aeqd +lat_0=90 +lon_0=0 to see a clear difference and create a layer called ‘iris_31_aeqd’.
Use the sf::st_crs function to guess the projection and sf::st_transform to change the projection.
#?st_crs
st_crs(iris_31)
par(mar = c(0,0,0,0), mfrow = c(1,2))
plot(st_geometry(iris_31))
iris_31_aeqd <- st_transform(iris_31, crs="+proj=aeqd +lat_0=90 +lon_0=0")
plot(st_geometry(iris_31_aeqd))Use map layers called ‘iris_31’ and ‘iris_31_aeqd’.
Units: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 51153.341 56509.22 53137.481 13756.98
[2,] 51153.34 0.000 20957.21 3086.011 49828.91
[3,] 56509.22 20957.212 0.00 11077.630 61863.18
[4,] 53137.48 3086.011 11077.63 0.000 54345.46
[5,] 13756.98 49828.910 61863.18 54345.458 0.00
Units: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 57206.922 62306.32 59385.205 13798.28
[2,] 57206.92 0.000 20957.48 3204.855 55368.23
[3,] 62306.32 20957.481 0.00 11076.220 66744.65
[4,] 59385.20 3204.855 11076.22 0.000 59843.30
[5,] 13798.28 55368.225 66744.65 59843.302 0.00
[1] FALSE
No, the two matrices are different.
The map layer called ‘iris_31’ contains the 5 digit codes of municipalities in its variable INSEE_COM and the 2014 population in its column P14_POP.
Use the classic functions of dplyr package: select, group_by et summarize. These functions also work with sf objects.
library(dplyr)
com_31 <- iris_31 %>%
select(INSEE_COM,P14_POP) %>%
group_by(INSEE_COM) %>%
summarize(P14_POP= sum(P14_POP)) %>%
st_cast("MULTIPOLYGON")
plot(st_geometry(com_31))The code of each municipality is not in the ‘sir_31’ database. To create it, you have to create a variable called INSEE_COM (5 digits) which concatenates the DEPET (2 digits) and COMET (3 digits) variables.
sir_31 <- readRDS("../data/sir_31.rds")
com_31 <- left_join(com_31,
sir_31 %>%
mutate(INSEE_COM=paste0(DEPET,COMET)) %>%
group_by(INSEE_COM) %>%
summarize(nb_of_rest= n()) %>%
st_set_geometry(NULL),
by=c("INSEE_COM"="INSEE_COM")
) %>%
mutate(nb_of_rest=ifelse(is.na(nb_of_rest),0,nb_of_rest))You have to use the database ‘table_MAUP.rds’ to have a match between the municipality code (CODGEO) and intercommunality code (EPCI).
table_MAUP <- readRDS("../data/table_MAUP.rds") %>%
select(CODGEO,EPCI)
epci_31 <- com_31 %>%
left_join(table_MAUP,by=c("INSEE_COM"="CODGEO")) %>%
group_by(EPCI) %>%
summarize(P14_POP=sum(P14_POP),nb_of_rest= sum(nb_of_rest)) %>%
st_cast("MULTIPOLYGON")
plot(st_geometry(epci_31))We can notice that the number of restaurants is very low or even equals to zero in most municipalities. We will therefore want to create a new map layer that partially uses the zoning of the municipalities and partially the EPCI zoning.
Start by creating the layers ‘EPCI_less100’ and ‘COM_more100’ corresponding respectively to the EPCIs which contain less than 100 restaurants and to the municipalities which belong to an EPCI containing more than 100 restaurants.
Use the do.call(rbind, list(EPCI_less100,COM_more100)) statement to merge these two layers.
EPCI_less100 <- epci_31 %>% filter(nb_of_rest < 100) %>%
setNames(c("territory","P14_POP","nb_of_rest","geometry"))
COM_more100 <- left_join(com_31,table_MAUP,by=c("INSEE_COM"="CODGEO")) %>%
filter(!EPCI%in%EPCI_less100$territory) %>%
select(-EPCI) %>%
setNames(c("territory","P14_POP","nb_of_rest","geometry"))
par(mar = c(0,0,0,0), mfrow = c(1,2))
plot(st_geometry(epci_31))
plot(st_geometry(EPCI_less100), col="lightblue",add=TRUE)
plot(st_geometry(epci_31))
plot(st_geometry(COM_more100), col="lightblue",add=TRUE)
mix_31 <- do.call(rbind, list(EPCI_less100,COM_more100))
plot(st_geometry(mix_31), col="lightblue")cartography package, simply add to each territory of this map a proportional circle layer related to the number of restaurants.The propSymbolsLayer function allows you to draw proportional circles.
library(cartography)
plot(st_geometry(mix_31), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
propSymbolsLayer(mix_31, var = "nb_of_rest", inches = 0.2)For the creation of ‘bks’ et ‘cols’, use the getBreaks et carto.pal functions of the cartography package. For the creation of the typo variable, you can use the cut function and apply the parameters digit.lab = 2 and include.lowest = TRUE.
library(sf)
library(cartography)
library(dplyr)
fra <- st_read("../data/fra.shp", quiet = TRUE)
mix_31 <- readRDS("../data/mix_31.rds")
mix_31 <- mix_31 %>% mutate(nb_rest_10000inhab = 10000 * nb_of_rest / P14_POP)
bks <- c(0,
getBreaks(v = mix_31$nb_rest_10000inhab[mix_31$nb_rest_10000inhab>0],
method = "quantile",
nclass = 3))
cols <- carto.pal("orange.pal", 4)
mix_31 <- mix_31 %>%
mutate(typo = cut(nb_rest_10000inhab,breaks = bks, dig.lab = 2, include.lowest = TRUE))cartography package, make the following map which contains in a choropleth layer the variable nb_rest_10000inhab and in a proportional circle layer the variable nb_of_rest. Do the same thing with the ggplot2 package.With cartography:
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(mix_31), border="azure")
plot(st_geometry(fra), col="ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)],add=TRUE)
choroLayer(mix_31, var = "nb_rest_10000inhab", breaks = bks, col = cols, border = "grey80", lwd = 0.4,
legend.pos = "topleft",
legend.title.txt = "Number of restaurants\nfor 10000 inhabitants",
add = TRUE)
propSymbolsLayer(mix_31, var="nb_of_rest", col=viridis::viridis(1,alpha=0.8),border=NA, legend.pos="left", legend.title.txt = "Number of restaurants", add = TRUE)
layoutLayer(title = "Restaurants", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)With ggplot2:
library(ggplot2)
map_ggplot <- ggplot() +
geom_sf(data = fra, colour = "ivory3",
fill = "ivory") +
geom_sf(data = mix_31, aes(fill = typo), colour = "grey80") +
scale_fill_manual(name = "Number of restaurants\nfor 10000 inhabitants",
values = cols, na.value = "#303030")+
geom_sf(data = mix_31 %>% st_centroid(),
aes(size= nb_of_rest), color = viridis::viridis(1,alpha=0.8), show.legend = 'point')+
scale_size(name = "Number of restaurants",
breaks = c(1, 500, 2000),
range = c(0,18))+
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(mix_31)[c(1,3)],
ylim = st_bbox(mix_31)[c(2,4)]
) +
theme_minimal() +
theme(panel.background = element_rect(fill = "azure",color=NA)) +
labs(
title = "Restaurants",
caption = "Insee, 2018\nKim & Tim, 2018"
)
plot(map_ggplot)mapview package. Try using different parameters to customize your map.For example, you can use the map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … parameters of the mapview function.
library(mapview)
library(sf)
library(cartography)
sir_31 <- readRDS("../data/sir_31.rds")
mapview(sir_31, map.types = "OpenStreetMap",
col.regions = "#940000",
label = paste(sir_31$L2_NORMALISEE, sir_31$NOMEN_LONG, sep = " - "),
color = "white", legend = TRUE, layer.name = "Restaurants in SIRENE",
homebutton = FALSE, lwd = 0.5)